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A new method has been developed to accelerate the convergence of explicit time- 
marching, laminar, Navier-Stokes codes through the combination of local preconditioning 
and multi-stage time marching optimization. Local preconditioning is a technique to 
modify the time-dependent equations so that all information moves or decays at nearly 
the same rate, thus relieving the stiffness for a system of equations. Multi-stage time 
marching can be optimized by modifying its coefficients to account for the presence of 
viscous terms, allowing larger time steps. We show it is possible to optimize the time 
marching scheme for a wide range of cell Reynolds numbers for the scalar advection- 
diffusion equation, and local preconditioning allows this optimization to be applied to the 
Navier-Stokes equations. Convergence acceleration of the new method is demonstrated 
through numerical experiments with circular advection and laminar boundary-layer flow 
over a flat plate. 


Nomenclature 

Roman letters 
A Cell area 

a , b Horizontal and vertical components of 
advection velocity, q 
c Speed of sound 

C Viscous flux Jacobian in the ^-direction 

C n Condition number, ratio of largest to 
smallest eigenvalues 

d Horizontal shift for Tchebyshev polynomial 

transformation 

E Energy per unit volume 

E Viscous flux Jacobian in the y-direction 

F Inviscid flux vector, ^-component 

G Inviscid flux vector, ^/-component 

H Total enthalpy per unit volume 

h Ratio of cell area to the length of the 
diagonal 

H Inviscid flux vector 

J Viscous flux vector 
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I Length from a cell corner to the diagonal 

along the advection velocity direction 
L 0 Reference length 

M Mach number 

p Pressure 

P Preconditioning matrix 

P n N th - order polynomial 

Pr Prandtl number, 0.72 

q Advection velocity 

Q o Reference velocity 

q XjV Heat flux components 

R Viscous flux vector, ^-component 

Rs Negative Real extent of the spatial 
operator’s Fourier Footprint 
Rt Negative Real extent of the temporal 
operator 

Re Reynolds number 

Re a Cell Reynolds number (generic) 

Re Ax Cell Reynolds number based on Ax 

Reh Cell Reynolds number based on h 

Re x Reynolds number based on ^-coordinate 
Res Residual 

S Conserved scalar quantity 

s Cell length 

S Viscous flux vector, y-component 

T Temperature 

T n iV th -order Tchebyshev polynomial 

U Blasius freestream velocity 
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u Velocity, horizontal component 

U Vector of conserved variables 

v Velocity, vertical component 

x Horizontal coordinate direction 

y Vertical coordinate direction 

z Complex number 

Subscripts 

( )o Reference quantity 

( ) E Euler 

( ), Cell face index 

( )j Cell index in y-direction 

( )ns Navier-Stokes 

( ) t Derivative with respect to time 

( ) Xj y Horizontal and vertical components or 
derivatives 

Conventions 

( ) Non-dimensional quantity (briefly) 

d( ) Differential quantity 

Symbols 

ak Runge-Kutta coefficient for stage k 

/3 Fourier frequency 

At Time step 

Ax Horizontal grid spacing 

Ay Vertical grid spacing 

£ Scaling for Tchebyshev polynomial 

transformation 

77 Blasius boundary layer coordinate 

7 Ratio of specific heats, 1.4 

5 Imaginary component 

k Parameter which controls upwind-biasing 

y Kinematic viscosity 

y , 0 Reference kinematic viscosity 

v Courant number 

< f) Flow angle with respect to the horizontal 

5ft Real component 

p Mass density 

g Prandtl-Glauert factor 

a Von Neumann number 

C Cell stretching parameter 

ov High frequency damping 

Tij Shear stress components 

6 Angle between h and the advection velocity 

£ High-frequency scaling for Preconditioner 

£ Transformed complex number 

Al Cell aspect ratio 

Ai q Equivalent cell aspect ratio in the flow 

direction 

T> Common denominator 


T Fourier Footprint 

Introduction 

T HIS work is motivated by the poor conver- 
gence of explicit time-stepping Computational 
Fluid Dynamics (CFD) codes for viscous flow prob- 
lems. The recent development of local preconditioning 
methods^' 3 ’ 1 and the flexibility of multistage time- 
marching schemes provide the tools to relieve this 
bottleneck.^ With the strong shift toward cache-based 
parallel architectures, implicit schemes are beginning 
to show limitations for this emerging environment. 

Typically an explicit time marching scheme is lim- 
ited by the minimum of an inviscid time step&S and a 
viscous time step . 0 The inviscid portion relates to the 
advective part of the equations and the viscous por- 
tion represents the dissipation component. In terms of 
Von Neumann stability analysis 0 these portions corre- 
spond to the extents of the Fourier Footprint^ along 
the imaginary axis and negative real axis, respectively. 

In the beginning of the “Euler era” (early 1980s), 
multistage time-marching schemes were optimized for 
the maximum imaginary extent of the stability do- 
main, given the number of stages of the scheme . 9 This 
allows inviscid problems to be solved efficiently on 
a single grid. Very soon, multigrid relaxation made 
its entrance,™ and the emphasis in multi-stage design 
shifted to choosing the coefficients such that maximum 
high-frequency damping results for a given number 
of stages, typically at half the maximum allowable 
time step.®®^ This optimization technique origi- 
nally was based on a scalar analysis, losing most of 
its validity when applied to a system of equations. 
The introduction of local preconditioning for the Euler 
equations,® which tends to equalize the advective time 
scales and concentrates the Fourier footprint, finally 
made it possible to optimize high-frequency damping 
for all types of waves admitted by the inviscid sys- 
tem.QEMS 

For a typical viscous flow problem, however, a ma- 
jority of the computational cells are limited in the size 
of their time step by the viscous criterion, due to the 
limited extent of the stability boundary along the neg- 
ative real axis. By making the multistage coefficients 
a function of the local cell Reynolds number, Re a, the 
stability boundary can be reshaped to alleviate this 
restriction .63’® 

For systems of equations, local preconditioning is 
necessary for equalizing various time scales admitted 
by the system.®®®®^ However, the analytical 
preconditioning analysis for the Navier-Stokes equa- 
tions is much more complicated than for the Euler 
equations, and is still far from complete. As a result, 

a Also, as a fringe benefit of local preconditioning, the solution 
accuracy at low Mach numbers is enhanced. 

b The locus of the Fourier transform in the complex plane. 
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the current paper uses a simple block-Jacobi method 
to extend the Euler preconditioner to viscous prob- 
lems . 1 ^ 

In the work of Lynn,®® multistage methods are 
optimized based on the precise Fourier footprint of the 
spatial operator. The optimization must be performed 
separately for each choice of cell Reynolds number. 
The multistage time steps can be tabulated as func- 
tions of Re a, or, more practically, given as functions 
that closely fit the data. Neither Lynn nor D. Lee 18 ac- 
tually carry out Navier-Stokes calculations with Re in- 
dependent multistage coefficients. 

In the present work the optimization procedure is 
reversed, and thereby greatly simplified. We start 
with a family of multistage schemes that has stabil- 
ity domains of desirable shape. Then, using the scalar 
advection-diffusion operator as a foothold, we fit the 
spatial footprint within the temporal stability bound- 
ary to achieve stability and possibly some other desir- 
able property, like prescribed high-frequency damping. 
Related work which optimizes the time step without 
consideration for the associated damping is given by 
Lorber et al E3E3 The transition from the advection- 
diffusion operator to the Navier-Stokes operator is 
accommodated by the local preconditioning. This pro- 
cedure provides a numerical relationship between Re a 
and the multistage time step coefficients. The result- 
ing ReA-dependent marching schemes are used to solve 
the advection-diffusion equation and compute laminar 
boundary-layer flow. 


equation 


St + a S x + b S y — —— ( S xx + S yy ) , (1) 

He 

where S is some conserved, scalar quantity; a and b are 
horizontal and vertical components of the advection 
velocity, respectively; and Re is the Reynolds number. 

Once we choose a discretization of the spatial terms, 
we are stuck with the resulting Fourier Footprint, so 
we will study it first. Afterward we will examine the 
temporal operator, and finally, the procedure used to 
fit one to the other. 

Spatial Operator 

For this study we will use central differencing for the 
viscous terms and the K-scheme®® for the advective 
terms. This yields, 


T = 




At 


Re Ax' 


(1 - cos (3 X ) 


(1 — k) aAt 
Ax 


(1 - cos 2f3 x ) 


4 

.aAt, . (k — 1) aAt . 

- t- 7 — sin + i — sin 2 p x 

Ax 4 Ax 


, , bAt „ At 

- (1 — k )— b 2 =■ 

' V Ay Re Ay 2 


(1 - cos P y ) 


(1 — k) bAt 
4 Ay 


(1 — cos 2 p y ) 


bAt . .(k — 1) bAt . 

— i — — sin p v + i — sin 2p v , 

Ay y 4 Ay y 


(2) 


Runge-Kutta Design 

Optimizing the time-marching scheme for efficiency 
or high-frequency damping is most easily accomplished 
by working in the complex space of the Fourier trans- 
form via Von Neumann stability analysis. 18 The over- 
riding goal, of course, is stability, i.e., to contain the 
Fourier Footprint of the spatial discretization inside 
the stability boundary of the time-marching scheme. 
A secondary goal is to provide some level of high- 
frequency damping. If only single-grid marching is 
desired, the tradition is to go for the largest time step 
possible. For multigrid marching it is necessary to 
maximize high-frequency damping, which means giv- 
ing up the maximum time step. The latter strategy 
actually turns out to also be more efficient for single- 
grid marching, as was shown by Tai. 113 Maximizing 
high-frequency damping without setting a target level 
for the damping, though, is not a good idea for a vis- 
cous equation solver, as this tends to reduce the time 
step to unacceptably low values.®® 

For the sake of simplicity we will start our analysis 
by examining the scalar model equation for the Navier- 
Stokes system of equations: the advection-diffusion 


where k £ [—1,1], providing a blending between fully 
upwind and central differencing, respectively. In a mo- 
ment we will look at the resulting Fourier Footprint, 
but first let us define some parameters, 


q 


\J a 2 + b 2 , 


0 

Al 


b 

= arctan - 


a 


Ax 
A y’ 


the advection speed, 
the flow angle, 

cell aspect ratio, 


and scale the footprint by fixing the position of the 
high frequency components (P x ,y = 7r) at — Rs- With 
these definitions, Eq. 0 becomes, 


— Rg = — 2A t 


(1 — k) — — (cos p + Al sin p) 
Ax 

,(! + Al 2 ) 


Re Ax 2 

In this form, two length scales have emerged, 

, Ax 

n = 


VI + Al 2 ’ 


c See Appendix [A| for non-dimensionalization. 


( 3 ) 

( 4 ) 
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the ratio of cell area to the length of the diagonal, and 


l = 


Ax 

cos <f> + Al sin (j> ’ 


( 5 ) 


which is the length from a cell corner to the diago- 
nal along the advection velocity direction. Figure |T] 
shows the cell geometry and the associated lengths. 
Moreover, h and l are related by l = h/ cos 9 where 9 


q 



is the angle between h and the advection velocity. 

Using the new length scale, h, we solve for At to fix 
the high frequency portion of the footprint at —R s , 


At = 


Rs 


2 

j[(l — k) cos 6 + 


2 

Reh 2 


( 6 ) 


Factoring out the q/h and converting the Reynolds 
number, Re , to the cell Reynolds number, Reh, we 
find, 


Rs h 

At = — -- 


Re h 


2 q Reh{ 1 ^ k) cos 9 + 2' 


( 7 ) 


Substituting this time step into Eq. 0 gives our scaled 
Fourier Footprint of the form, 


T = - 


Rs h 


Reh 


2 Ax Reh cos 9 
Re Ax cos <j) + 2 


Reh 


(1 — cos f3 x ) + i cos <j> sin f3 x 


M 


f Re Ax sin cf> + 2 




Reh 


(8) 


(1 - COS Py) 

iAl sin 4> sin (3 y 


Finally, we reveal the Fourier Footprint in Fig. 0 which 
shows the effect of various k’s for Rs = 1, tR = 1, 
and Re = 1. Defining a non-dimensional time step 
for the advection-diffusion equation is not trivial, see 
Appendix [B] for more details. 

Now that we have the spatial operator’s footprint, 
we turn our attention to the temporal stability bound- 
ary. 
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a) Upwind (k. = —1). 
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c) Third-order in ID only d) Central difference 


(k=1/3). 


(k, = 1). 


Fig. 2 Effect of n on the footprint of the higher-order 
scheme (yft = 1, <t> = 45 deg, Re h = 1, /3 x ,f3 y £ [0,27r], 
A/3 = 2tt/40). 


Temporal Scheme 

For this study the stability boundary is given by a 
transformation of Tchebyshev polynomials attributed 
to Manteuffel, 25 


Pn(z) 


T n (^) 
Tn(D ’ 


( 9 ) 


where n is the order of the respective polynomials, z is 
a number in the left complex plane, and the sequence 
of Tchebyshev polynomials are given by recursion, 


r 0 (C) = i, (io) 

Ti(0 = C, 

T„ +1 (C) = 2 C T n ( C) - T„_a(C), n > 1. 

The resulting coefficients of P n (z) are chosen such that 
the \P n (z) | = 1 stability boundary remains simply con- 
nected and the blended polynomial is scaled properly, 
i.e., 

Pn{ 0 ) = 1 , ( 11 ) 


and 


dP n 

dz 


= 1. 


(12) 


z=0 


For example, a four stage scheme would take the 
form, 


. , 32d 3 — 16e 2 d 

Pa\z) = 1 — z 


V 

48 d 2 - 8s 2 
V 


o 32 d, o 8 4 
~ - + —z 4 , 

V V ’ 


(13) 


where V = 8 d 4 — 8 d 2 e 2 + e 4 . To satisfy the conditions 
given by Eqs. 0 and 0 we find, 

£= y / 4d(2+d) + v / 16d 2 (2+d) 2 -8d 3 (dTiy. ( 14 ) 
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Thus, we now have a one-parameter family of stability 
boundaries governed by d. 

Next, we recall that a 4-stage Runge-Kutta scheme 
has the form of 

Pa{z) = 1 + 04 z + a 4CK3 z 2 + 040:30:2 £ 3 + 04030201 z 4 . 

(15) 

So by equating the coefficients in Eqs. [13] and |T0 we can 
determine the coefficients, «&, for the corresponding 
Runge-Kutta scheme. 

Fitting 

The procedure employed is analogous to blowing up 
a balloon into a bear trap: a given stability boundary 
(the bear trap) is chosen; next, a given high frequency 
damping level is located on the negative real axis; 
and then the cell Reynolds number is increased until 
some part of the spatial operator’s Fourier Footprint 
(the balloon) reaches a specified damping contour of 
the temporal operator. Figure 0 depicts the parame- 
ters, the stability boundary, and the spatial operator 
Fourier Footprint. 



Fig. 3 Sketch of optimization procedure. 


The fitting algorithm is as follows: 

1. Pick a value for d. This sets the shape and extent 
of the Tchebyshev stability boundary. 

2. Compute — Rt , the negative Real extent of the 
Temporal operator. (Hint: Rt = 2 d.) 

3. Given the negative Real extent of the stability do- 
main, — Rt , as a starting point, find the negative 
Real extent of the spatial operator, —R s , which 
satisfies the prescribed high frequency damping ,0 
ov, by moving toward the origin along the Real 
axis, computing |P„(,s)|. 

4. Using the viscous limit ( Reh — > 0) as an initial 
guess increase the cell Reynolds number until the 
prescribed damping conditions for the entire foot- 
print or its high-frequency part are met. 

d Or, if unable to find the requested damping, use the maxi- 
mum damping available. 


Figure 0 on the following page shows a sequence of 
stability plots for different values of d. Superimposed 
are the spatial operator Fourier Footprints which yield 
the maximum time step given a high-high frequency 
damping of 0.5. The contours in Fig. 0 represent levels 
of the temporal operator’s amplification (or damping) 
factor, |P|, from 0.1 to 1.0 in increments of 0.1. 

Notice the large negative extent of the stability do- 
main along the real axis, made possible by the Tcheby- 
shev polynomial. The domain maximally extends to 
— 2n 2 for an n-stage scheme; however, for this choice 
the scheme is only stable in the limit of Re a = 0 (refer 
to Fig. |'l(d)| ). 

Figure 0 on page [7] shows the reference time step 
and Runge-Kutta coefficients as a function of the 
cell Reynolds number for the resulting 4-stage Runge- 
Kutta scheme. This provides the link between the cell 
Reynolds number, the time step, and the multistage 
coefficients for a given cell. Note that the time step 
maintains an appreciable value all the way down to 
Reh = 0.1, which is smaller than many people use to 
resolve the boundary layer in a Navier-Stokes calcula- 
tion. 

Extension to Systems (Local Preconditioning) 

Local preconditioning is a technique to remove 
stiffness from a system of equations. In the 
aeronautical CFD community, the set of time- 
dependent, Reynolds-averaged, compressible, Navier- 
Stokes equations — sometimes simplified to the (invis- 
cid) Euler equations — are typically solved in an iter- 
ative fashion to achieve a steady-state solution. Al- 
though in most cases only the steady-state solution 
is desired, the time-dependent equations are still em- 
ployed so that the marching problem is well posed 
for all Mach numbers. The time-dependent Euler 
equations are hyperbolic and admit real wave-like fun- 
damental solutions; the Navier-Stokes equations are 
mixed parabolic-hyperbolic and admit damped travel- 
ing waves as well as non-propagating damped modes. 

Convergence to steady-state for inviscid calculations 
is impaired in some flow regimes due to the spread in 
the characteristic wave speeds. The spread is largest 
for waves moving in the streamwise direction, in which 
case their speeds are q , q + c, and q — c, where q is the 
total flow speed and c is the speed of sound. The 
ratio of the largest to smallest wave speed is termed 
the condition number, C n , and serves as a measure of 
“stiffness”. Figure 0 on page [7] shows the condition 
number as a function of Mach number for the Euler 
equations. The solid line represents the original Euler 
equations and it is apparent that infinite stiffness oc- 
curs in both the subsonic and transonic regimes; as the 
Mach number further increases, the condition number 
asymptotes to the ideal value of one. 

The dashed line in Fig. 0 represents the Euler equa- 
tions preconditioned with the Van Leer-Lee-Roe pre- 
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5R 


e) d = —16 

Fig. 4 Optimized stability plots as a function of the temporal operator parameter, d, for Fromm’s Scheme (k. = 0) 
with a prescribed high frequency damping of 0.5 with contours of damping every 0.1 and the size of the frequency 
symbols scaled by their respective damping. 


6 of D 


American Institute of Aeronautics and Astronautics Paper 99-3267 


10 


(16) 



Fig. 5 The reference time step and Runge-Kutta co- 
efficients as a function of the cell Reynolds number for 
the fitting presented in Fig. g|. 



M 


Fig. 6 Condition number as a function of Mach number 
for the Euler equations. 


conditioner.® This preconditioner comes closest to 
achieving equalization of wave speeds, without affect- 
ing their direction of propagation or their hyperbolic- 
ity. The dashed line shows it is possible to eliminate 
completely the stiffness as the Mach number goes to 
zero, greatly reduce the transonic stiffness, and, in gen- 
eral, substantially lower the condition number for the 
system of equations. This makes the system behave as 
a scalar equation which has only a single wave speed, 
so we can apply the Runge-Kutta scheme developed in 
the previous section after accounting for some scaling 
factors. 

The inviscid portion of local preconditioning used 
in this paper is that of Van Leer-Lee-Roe 01 with a 
modification according to D. Lee^9’E3 and Lynn. 10 
In stream-aligned, hyperbolic symmetrizing variables, 


[dp/ (pc), dn, dv, dp — c 2 dp ] T , it has the form, 


r f m 2 

? g2 



0 


0 O' 
1+Jr 0 0 
0 0 
0 0 1 . 


with g = sj |1 — M 2 1 and £ = g/(g + At,), where M. q 
is an effective streamwise aspect ratio given by sum- 
ming the appropriate velocity projections for each cell 
face, i, 


= J2 l \ uAx i + vA Vi\ 

9 J2i\ vAx i - uA Vi\' 


(17) 


This form is due to Darmofal. 123 

Preconditioning of the Navier-Stokes equations is 
accomplished according to the method proposed in 
Ref. |2I, namely, the addition of the viscous Jacobians: 

p ~‘s = p » ,+ Mb) |C+ “»’ (18 > 


where C and E are the Jacobians of the viscous fluxes 
in the x and y directions, respectively. 

Thus, instead of marching to the steady state with 


U t =Res( U), (19) 

march with a preconditioned residual 

U t = P(U)i?es(U). (20) 


This can be viewed as marching with a matrix time 
step as opposed to a scalar time step. However, the 
following three elements of an existing code need to 
be modified to account for the new, preconditioned 
system: 

• Time-step definition 

• Artificial dissipation 


• Boundary conditions 

First, since preconditioning serves to collapse all 
wave speeds to the total flow speed, q, the time step 
should be based on the flow speed, q , and not the 
largest acoustic wave speed, q + c, which is typically 
used in the unconditioned case. 

Next, for an upwind scheme the artificial dissipation 
matrices |A| and |B| need to be evaluated in terms of 
preconditioned quantities, i.e., they become P -1 |PA| 
and P _1 |PB| for the preconditioned system. 

And lastly, if one is employing “weak” (image cell) 
boundary conditions and using the modified upwind 
scheme to solve for the boundary fluxes, no change 
is necessary for the boundary conditions. However, if 
one is using explicit characteristic boundary conditions 
and/or employs “strong”, specified boundary fluxes, 
these procedures need to be modified to account for the 
preconditioned equations that are now being solved. 
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Numerical Experiments 

Two numerical experiments were performed, the 
first on the scalar model equation for the Navier-Stokes 
equations: advection-diffusion; and the second experi- 
ment solved the Navier-Stokes system of equations for 
boundary layer flow over a flat plate. 

Circular Advection 

To remove the added complication of precondi- 
tioning a system of equations, the scalar advection- 
diffusion model problem was first considered. The 
case considered here is that of circular advection in 
the upper half plane, i.e., the domain consisted of 40 
equally-spaced cells along the x direction from -1 to 1 
and 20 equally-spaced cells along the y direction from 
0 to 1, with a top-hat profile entering at the lower left 
rotating clockwise to exit at the lower right. Figure 0 
on the next page shows a representative solution for a 
Reynolds number of 500. 

The number of iterations to converge the L 2 er- 
ror norm to 10 ~ 10 are recorded in Table 0 for a 
host of Reynolds numbers for both the variable- 
coefficient Runge-Kutta scheme and the constant- 
coefficient scheme. The fixed-coefficient scheme uses 

Table 1 Comparison of constant-coefficient Runge- 
Kutta scheme to variable-coefficient Runge-Kutta 
scheme for Circular Advection on a 40 X 20 grid. 


Iterations 

Re 

Variable Fixed 

1 

868 

4,271 

10 

422 

2,736 

25 

252 

1,374 

50 

151 

736 

100 

110 

387 

500 

85 

136 

1,000 

77 

117 

5,000 

81 

98 

10,000 

82 

99 


Table 0 clearly demonstrates the devastating ef- 
fects that low Reynolds numbers (and hence low cell 
Reynolds numbers) have on the standard advection- 
optimized time marching scheme. As anticipated by 
the Fourier analysis, the variable-coefficient scheme re- 
tains a healthy convergence rate even with local cell 
Reynolds numbers well below 0.01. 

Laminar Boundary Layer 

The flow solver for this portion of the study was pur- 
posely kept as simple as possible to isolate the effect 
of the variable Runge-Kutta coefficients. The solver 
consists of a cell-centered, finite- volume-based scheme 
which uses Roe’s Flux Difference Splitting^ for the 
inviscicl fluxes and central differencing for the viscous 
terms on a structured grid of quadrilaterals. In inte- 
gral form, 

/ U t dA+ ® H ds = <J J ds, 

J A Js Js 

where U is the conserved state vector defined as 
(p, pu , pv, pE) T , H is the inviscid flux vector, 

(F? + G j) ■ ds, 

and J is the viscous flux vector, 

(Ri + Sj) • ds. 

The Cartesian components of the inviscid flux are 
F = (pu, pu 2 + p, puv, puH ) T , 

and 

G = (pv, pvu, pv 2 + p, pvH) T , 
while the viscous components are 

R (0, T'xx, X x y, VT X x~\~VT X y q x ( , 

and 

S (0 , Xxy , Tyy , UT X y ~\~VTy y (jy ) 


advection-optimized Runge-Kutta coefficients^ due to 
Tai 113 for maximizing damping over frequencies from 
7r/4 to 7r while the variable-coefficient scheme uses the 
cell Reynolds number-dependent coefficients that were 
developed in the preceding section. Note that for this 
problem the average cell Reynolds numbers are a fac- 
tor of 20 less than the Reynolds number and actually 
tend toward zero near the center of the advection ve- 
locity field since Ren = Re h q. However, for very low 
values of q we found that the local time-step given by 
Eq. 0 on page 0 must be capped to avoid numerical 
instability regardless of the marching method used. 

Coefficients: ai=0.2131, (*2=0.4364, 03 = 0. 7641 with 

u=1.1727. 


The governing equations are non-dimensionalized by 
freestream speed of sound, c, and density, p, so that 
the viscous stresses and heat flux terms are given by 

M 2 

Tex — ■^ x ^2/' 7 

M 2 

T yy ~ Re ^3 ^ Vy ~ Ux ’ 

M , 

Try n P (lly 4” ) > 

He 


M 


0.x 

Qy 


(7 — l)i?ePr 
M 

(7 — l)i?ePr 


pT x , and 


pT y . 
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Fig. 7 Circular advection for a Reynolds number of 500 including velocity vectors and solution contours every 
0.1 from 0.05 to 0.95. 


The system is closed with an equation of state, 


P = (7 - 1) 


E- 


+ v 2 * ) 


where E is the total energy per unit volume, H is the 
total enthalpy per unit volume, T is the temperature, 7 
is the ratio of specific heats (1.4), /r is the viscosity, Pr 
is the Prandtl number (0.72), and Re is the Reynolds 
number. 

The numerical test problem chosen is that of a lam- 
inar boundary layer flow; specifically, two-dimensional 
subsonic flow over a flat plate. The unit-length 
Reynolds number is 10,000 and the freestream Mach 
number is varied between 0.05 and 0.3. The compu- 
tational domain and mesh are shown in Fig. || on the 
following page. The plate is 4 units long, with the up- 
stream boundary 2 units away from the leading edge 
and the upper boundary is placed at 1.2 units. There 
are 36 cells evenly distributed along the cc-direction 
(24 cells on the plate and 12 upstream) and 40 cells 
exponentially stretched in the y-direction according to 

eM - 1 - 1 . i 

Vj = y nj ^ _ 1 for j = l,... no- 
where the wall spacing is set to have an acoustic 
cell Reynolds number of 10, yielding a stretching pa- 
rameter, c, of 1.1360.0 By employing characteristic 
boundary conditions, neither second-order boundary 

'Parabolic stretching, as used in Refs. £jo| and [ijj, was origi- 
nally tried until, upon further examination, it became apparent 
that it had a canonical first-point stretching factor of three! 


conditions®’® nor solution-assumed grids®’®’® are 
necessary to capture the boundary layer gradients (see 
Fig. 0 on the next page where the computed results^ 
are compared to the Blasius results^]). 

Table 0 shows the number of iterations to converge 
the L 2 error norm to 10 ~ 6 for a range of Mach num- 
bers for both the fixed-coefficient Runge-Kutta scheme 
and the variable-coefficient Runge-Kutta scheme.0 The 

Table 2 Comparison of constant-coefficient Runge- 
Kutta scheme to the variable-coefficient Runge-Kutta 
scheme for boundary layer flow over a fiat plate. 


M 

Fixed 

Uncond. Precond. 

Variable 

0.3 

21,322 

7,667 

3,102 

0.1 

45,371 

7,482 

1,945 

0.05 

97,263 

8,231 

1,543 


fixed-coefficient scheme was run with and without pre- 
conditioning. The results clearly show the benefits of 
preconditioning, with local preconditioning providing 
nearly Mach number independent convergence. How- 
ever, we note some degradation as the cell Reynolds 
numbers become lower as the Mach number is lowered. 
This effect was predicted by D. Lee®’® and shows 

g The unconditioned and preconditioned solutions are indis- 
tinguishable. 

h Strictly speaking, Blasius is only valid for incompressible 
flow, but up to Mach 0.3 compressibility effects are negligible. 

1 CPU timing data are not presented for this proof-of-concept 

study since no optimization of the new implementation has been 

attempted. 
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Fig. 8 Computational domain and grid. 



a) u-velocity b) u-velocity 

Fig. 9 Boundary layer solution profiles at X/L = 0.521 for Mach 0.1 flow. 


the limitation of simply using the viscous Jacobians 
to augment the Euler preconditioning (see Eq. |lS| on 
page 0 ). 

Most important, the results for the variable- 
coefficient scheme show even further improvement in 
convergence efficiency: the variable-coefficient scheme 
reaches convergence with nearly a factor of five fewer 
iterations than the fixed-coefficient scheme. 

Conclusions 

A four-stage time-marching scheme was optimized 
for a viscous model problem. The Runge-Kutta coef- 
ficients as a function of cell Reynolds number were 
designed to yield a specified damping level while 
maintaining large time steps. These variable coeffi- 
cients were applied to a two-dimensional, cell-centered, 
Navier-Stokes solver which incorporated local precon- 
ditioning. The results for a subsonic flat plate bound- 
ary layer flow indicate that by employing local Runge- 
Kutta coefficients as a function of cell Reynolds num- 
ber, convergence to steady state is greatly enhanced. 

Since this is a proof-of-concept study, many other 
variables were purposely held constant so that the 


effects of the cell Reynolds number-dependent Runge- 
Kutta coefficients could be ascertained. Given the 
qualitative improvement over the standard method, 
the new approach appears to justify further develop- 
ment. 
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Appendix 

A Non-dimensionalization of the 
Advection-Diffusion Equation 

We begin with 

St + Cl S X + b Sy = fl(S XX + Syy ) , (A.l) 
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then define non-dimensional quantities, 


5 S - a z b - V 

S = a = 7 V’ b=z 7V’ 

Wa Oo no 

i= b s= b mdi= i^- 

Substituting into Eq. |A.1| yields, 

d(S 0 S ) d(S 0 S) ^{SpS) 

+ Qoa-^—r + Qob 


d(L 0 /Q 0 t) 

= MoM 


(A.2) 


d{L 0 x) d(L 0 y) 

( d 2 (S 0 S) d 2 (S 0 S) \ 

\8LpxdLpX dLpydLpy) ' 

Collecting constant terms gives, 


(A.3) 


SpQo^ . S 0 Qp_- S 0 Q o-- 

— — bi H - — ab x H — bby 

Lp L 0 Lp 

= 7n ^° fi (, S ^ + Syy ) , (A. 4) 

L o 

while multiplying through by Lp/(SpQp) and defining 
Re = SpLp/Ho yields, 

Sf + aSx + bSy = (S xx + Syy ) . (A.5) 

If we chose Ho = H and drop the bar notation, we have 
simply, 

St + a S X + b Sy = J^{S XX + Syy ) . (^-6) 

B Time-Step Normalization for the 
Advection-Diffusion Equation 

Let us try to find a non-dimensional time step that 
includes both advective and diffusive effects. 


where k £ [—1,1]. 

Typically one defines parameters according to the 
advective and diffusive terms, e.g., 


v 


aAt 

Ax 


and 


At 

~ Re Ax 2 ’ 


where v is the Courant number 13 and a is the Von 
Neumann number.™ This gives a footprint in Fourier 
space as shown in Fig. |B.1| which is an ellipse defined 
by the Courant number and the Von Neumann num- 
ber. (The higher-order scheme yields an ellipse only 
for k = 1; in general it is egg-shaped.) 





Fig. B.l Fourier footprint for the one-dimensional 
Advection-Diffusion equation using first-order upwind 
spatial discretization for the advective term and central 
differencing for the diffusion term (/3 £ [0,27 t]). 


Now, suppose we want to fix the location of the neg- 
ative real extent of this footprint. We can do this by 
examining the high frequency limit (/? = 7 r) of the 
Fourier footprint given in Eq. |l > . 2| . 


-Rs = -2 


/ aAt 

\ A:r 


2 


At \ 
Re Ax 2 ) ’ 


(B.4) 


where —Rs is the specified negative Real axis extent. 
Now we simply solve for the time step, At, 


One-Dimension 

To provide a simple foundation, we begin with the 
one-dimensional case, 

S t + aS x = ±-S xx , (B.l) 

Re 

where S is some conserved scalar quantity, a is the 
advection speed, and Re is Reynolds number. Cen- 
tral differencing for the diffusive term and first-order 
upwind discretization of the advective term yields, 


T=- 


( aAt 
\ Ax 


2 


At \ 
Re Ax 2 ) 


(1 — cos /3) — i sin f3 , 
(B.2) 


where (3 represents the Fourier frequency. Alterna- 
tively, employing a higher-order «;-scheme^’E3 for the 
advective term yields, 


T = 




At 


(1 — cos/3) 


— 1- 


(k — 1) aAt 
4 Ax 
. (3 — k) aAt 


Re Ax 2 , 

(1 — cos 2/3) (B-3) 

(1 — k) aAt 


Ax 


sin/3 - 


Ax 


sin 2/3 , 


At = 


Rs 

2 

a I 2 
Ax ' Re Ax 2 


For the K-scheme the result is similar, 


At = 


As 

2 


(1 - k) 


a 

Ax 


2 

Re Ax 2 


(B.5) 


(B.6) 


Thus, by specifying such a time step, the high fre- 
quency components of the Fourier footprint remained 
fixed to —Rs- We could also write this in terms of 
a cell Reynolds number, 0 Re a x = (clAx)/h where 
H = 1 /Re as shown in Appendix |A]. For the first-order 
discretization of the advective term this yields, 


or 


At = 


a 

Ax 


Rs_ 

2 



At - — — ReAx 
2 a Re^ x , + 2 


while the K-scheme gives, 

_ R s Ax Re Ax 

2 a Re& x {l — re) + 2 


(B.7) 


(B.8) 


(B.9) 
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Alternatively, we could express these as conditions on 
the Courant and von Neumann numbers, 


Rs 


2(v + 2a) first-order 

2[(1 — k)v + 2a] higher-order 


(B.10) 


Courant numbers and Von Neumann numbers, e.g., 
aAt At 


Ax ’ 

bAt 
Ay ’ 


and a„ = 


Re Ax 2 ’ 
At 


ReAy 


2 ■ 


Now, we insert our definition of At from Eq. |B.7| into 
Eq. |B.2| on the page before to verify that the footprint 
is now only dependent on Rg and the cell Reynolds 
number, 


T=- 


Rs 

2 


(1 — cos /?) + i 


Re Ax 
Re Ax + 2 


sin /3 . 


(B.ll) 


So Rs controls the negative real extent and the cell 
Reynolds number dictates the ellipticity. 

Two-Dimensions 

While the one-dimensional case is relatively 
straightforward in terms of parameter choices, etc., 
adding another dimension creates some “pseudo” am- 
biguities which are not typically resolved in the CFD 
community. 

The governing equation is now, 

S t + aS x + bSy = ^ ( S xx + S yv ). (B. 12) 

Again, using a first-order, upwind approximation for 
the advective terms and central differencing for the 
diffusion terms yields, 


_ ( aAt At 

F = — l — b 2 — ^ 

V Ax Re Ax 2 

/ bAt o At 
~ V Ay + 2 ReAy 2 


aAt 

(1 - cos /3 X ) - i sin fj x 

(1 — cos j3 y ) — sin Py ■ 
(B.13) 


Now there are two Fourier frequency components, [3 X 
and (3 y . The corresponding higher-order result is, 


_ . . .aAt At 

T = - (1 — K) — b 2 


Ax Re Ax 2 

(1 — re ) aAt 


(1 - cos/3*) 


4 Ax 


(re — 1) aAt . 

+ — : t— sm 2/3* 

4 Ax 

/C , bAt „ At 

- (1 — k) —r b 2 


V Ay 1 “ ReAy 2 

(1 — k) bAt . 

4 A y ^ ~ C ° S 


(1 - cos (3y) 


(B.14) 


bAt . n (re — 1) bAt . „ „ 

— i — — sm Ijy + i — sin2A. 

Ay v 4 Ay y 


Instead, we propose to express the flow quantities in 
cylindrical coordinates and the geometric quantities in 
terms of Ax and the cell aspect ratio, 

q = \J a 2 + b 2 , 4> = arctan -, 

a 


and Al = — — . 

Ay 

( Cf . Fig. |T] on page f| in the main text.) Applying 
these definitions to Eq. gives, 


T=- 


qcoschAt „ At , 

A l-cos^ x 

Ax ReAx 1 

. q cos (f>At 


Ax 


sin f3 x 


_ + 2 ^ M 


Ax 

_.^^sin0At sto 

Ax 


w (1 " cos/3 » ) 


(B.15) 


As was done for the one-dimensional case, we ex- 
amine the high frequency limit (/3 X = (3 y = tt) of the 
Fourier footprint given in Eq. |B.15| . 


— Rs = — 2At 


cos 4> + AH sin </> 
9 Ax 


2 1 + Al 2 \ 
Re. Ax 2 ) ; 

(B.16) 


and, solving for the time step such that the nega- 
tive Real extent of the footprint remains fixed at —Rs 
yields, 

Rs 


At = 


2 

q cos 6 2 ’ 

h ' Reh 2 


(B.17) 


where h is a new length scale that corresponds to the 
ratio of cell area to the cell’s diagonal length and the 
angle 9 is the angle between h and the flow direction. 
(Again, cf. Fig. |T] on page |j.) The re-scheme yields a 
similar form, 


At = 

(1 - re) 


Rs 

2 

q cos 6 | 
h h 


2 

Reh 2 


(B.18) 


Again, as was done for the one-dimensional case, we 
can re-write this in terms of a cell Reynolds number, 
Reh = ( qh)Re , so Eq. becomes, 


From here there are many avenues to follow, one 
of the most prevalent is to define component-wise 


Rs h R e h 

2 q Reh cos 9 - b 2 


(B.19) 
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and the K-scheme yields a similar result, 


^ _ Rs_ h Reh 

2 q Reh{ 1 — k) cos 0 + 2 ' 


(B.20) 


Note the differences with respect to Eq. |B.8| and 
Eq. |B.U| on page 0 In the two-dimensional case we 
have a new length scale, h, which amounts to some- 
thing akin to the harmonic average of the two length 
scales, 



(B.21) 


Ay 2 


This favors the smaller of the two components, h can 
also be expressed as the ratio of the cell area to the 
length of the diagonal, 


Ax Ay 

\J Ax 2 + Ay 2 


(B.22) 


Another twist is that the flow speed component in the 
“/i” direction is what governs the advective portion 
and not the full flow velocity as one might anticipate. 

Now we will substitute our value of At into the orig- 
inal Fourier footprint of Eq. mm on the preceding 
page, 


T= - 


Rs h 


Re h 


2 Ax Reh cos 6 
Re Ax cos (j)+ 2 


Ren 


M 


(1 — cos fi x ) + i cos ^ sin fi. 
f Re Ax sin <j) + 2 


V 


Reh 

+ iA\ sin cj) sin fi y 


(1 — COS fiy ) 

(B.23) 


This form is not as straightforward as the one- 
dimensional results, but examining the case of Al = 1, 
<f> = 45 deg, we do recover a similar form; 


J _ (cos fix + COS fiy ) 

2 

. Re h (sin fix + sin fiy ) 

+ i — 

Re h + 2 2 



(B.24) 


which for fi = fi x = fi y , gives the identical result to 
Eq. |B.11| on the page before. Figure @ on page [| in the 
main text shows Fourier footprints for four values of n 
for the higher-order scheme with Al = 1, <j> = 45 deg, 
and Reh = 1. 
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